%% Initialization
%clear all
close all
clear global
clc
delete(gcp('nocreate'))
distcomp.feature( 'LocalUseMpiexec', false )

% restoredefaultpath;
rehash toolboxcache;

MyPool = parpool;
%% Add paths
addpath('../Matlab')

%% Read network (GRN+interface nodes)
G = importdata(strcat('adjmat-contextualized-',subfix,'-interface.txt'));
InData = importdata(strcat('pheno-contextualized-',subfix,'-interface.txt'));

%Formatting
RowName = G.textdata;
RowName(1) = [];
G = G.data;
BoolData = InData.data;
GNames = InData.textdata;

initial_state = 1; % from experimental data

% Nodes indexes
delim = '\t';
to_perturb = importdata(strcat('to-perturb-',subfix,'.txt'),delim);
node_ind = zeros(size(to_perturb));
for i = 1:size(to_perturb)
    node_ind(i) = find(strcmp(GNames, to_perturb(i)));
end

% State of nodes to perturb
pert_BoolData = BoolData(node_ind, initial_state);
ups = node_ind(pert_BoolData==1); %index in BoolData(pert_ind, initial_state)
downs = node_ind(pert_BoolData==0);
both = node_ind(pert_BoolData==2);

%% Unperturbed attractor
% Unperturbed network (no nodes connected)
GRN_ind = ~ismember(1:size(G,1), both); %logical vector to select only GRN
unp_G = G(GRN_ind, GRN_ind);
unp_BoolData = BoolData(GRN_ind,:);

% Unperturbed attractor
[Att_states, count] = boolSim_pert(unp_G,unp_BoolData(:, initial_state),1000,[]);
m = [sum(Att_states==unp_BoolData(:, initial_state)) / size(Att_states, 1)*100, ... 
       sum(Att_states==unp_BoolData(:, initial_state)), count, sum(~(Att_states == unp_BoolData(:, initial_state)))];
unperturbed_Att_states = Att_states;
m_flipped = strjoin(GNames(find(Att_states ~= unp_BoolData(:, initial_state))),';');

%% Perturbations
results=[];
perturbations = char.empty(0);
M_flipped = char.empty(0);

parfor i=1:(size(node_ind,1)-1)
    i_n = node_ind(i,:);
    add_i = ismember(i_n,both); % do I need to add it, or is it already in the GRN, and with a fixed state? 1=add,0=fixed
    
    for j= (i+1): size(node_ind,1)
        j_n = node_ind(j,:);
        add_j = ismember(j_n,both); % do I need to add it, or is it already in the GRN, and with a fixed state?
        
        for w = (j+1):size(node_ind,1)
            w_n = node_ind(w,:);
            add_w = ismember(w_n,both); % do I need to add it, or is it already in the GRN, and with a fixed state?
            
            for z = (w+1):size(node_ind,1)
            z_n = node_ind(z,:);
            add_z = ismember(z_n,both); % do I need to add it, or is it already in the GRN, and with a fixed state?
        
                if (~add_j && ~add_i && ~add_w && ~add_z) % if all are already in the network
                    ind_pert = [i_n j_n w_n z_n];
                    [Att_states, count] = boolSim_pert(unp_G,unp_BoolData(:, initial_state),1000,ind_pert); 
                    %Compare with unperturbed attractor
                    expected_Att_states = unperturbed_Att_states;
                    expected_Att_states(ind_pert) = ~expected_Att_states(ind_pert);
                    flipped_genes = sum(~(Att_states == expected_Att_states));

                    temp_res = [sum(Att_states==expected_Att_states) / size(Att_states, 1)*100, ... 
                        sum(Att_states==expected_Att_states), count, flipped_genes]; 
                    temp_flipped = strjoin(GNames(find(Att_states ~= expected_Att_states)),';');
		    % save name of perturbation
                    j_name = strcat(GNames(j_n),'_',num2str(1*~unp_BoolData(j_n,initial_state)));
                    i_name = strcat(GNames(i_n),'_',num2str(1*~unp_BoolData(i_n,initial_state)));  
                    w_name = strcat(GNames(w_n),'_',num2str(1*~unp_BoolData(w_n,initial_state)));
                    z_name = strcat(GNames(z_n),'_',num2str(1*~unp_BoolData(z_n,initial_state)));
                    temp_name = char(strjoin({i_name{:},j_name{:},w_name{:},z_name{:}},';'));

                elseif (add_j && add_i && add_w && add_z) % if I have to add all
                    ind_pert = [i_n j_n w_n z_n];
                    pert_G = G(~ismember(1:size(G,1), both(~ismember(both,ind_pert))), ~ismember(1:size(G,1), both(~ismember(both,ind_pert))));
                    pert_BoolData = BoolData(~ismember(1:size(G,1), both(~ismember(both,ind_pert))),:);

                    n=4; % perturbed nodes
                    k = dec2bin(0:(2.^n)-1)-'0'; % possible combinations

                    temp_res = zeros(size(k,1),4);
                    temp_name = cell(size(k,1),1);
		    temp_flipped = cell(size(k,1),1);
                    for l=1:size(k,1)
                        p_BoolData = pert_BoolData;
                        p_BoolData(p_BoolData(:,initial_state) ==2,initial_state) = k(l,:); % assign combinations of values
                        % Simulation
                        [Att_states, count] =boolSim_pert(pert_G,p_BoolData(:, initial_state),1000,[]);
                        %Compare with unperturbed attractor
                        comparable_Att_states = Att_states(GRN_ind,:);
                        flipped_genes = sum(~(comparable_Att_states == unperturbed_Att_states));
                        temp_res(l,:) = [sum(comparable_Att_states==unperturbed_Att_states) / size(Att_states, 1)*100, ... 
                            sum(comparable_Att_states==unperturbed_Att_states), count, flipped_genes]; 
			temp_flipped{l,:} = strjoin(GNames(find(comparable_Att_states ~= unperturbed_Att_states)),';');

                        i_name = strcat(GNames(i_n),'_',num2str(k(l,1)));
                        j_name = strcat(GNames(j_n),'_',num2str(k(l,2)));
                        w_name = strcat(GNames(w_n),'_',num2str(k(l,3)));
                        z_name = strcat(GNames(z_n),'_',num2str(k(l,4)));
                        temp_name{l,:} = char(strjoin({i_name{:},j_name{:},w_name{:},z_name{:}},';'));
                    end          

                else % if I have to add only some
                    indeces = [i_n j_n w_n z_n];
                    ind_interf = indeces([add_i,add_j,add_w,add_z]);
                    ind_pert = indeces(~[add_i,add_j,add_w,add_z]);

                    pert_G = G(~ismember(1:size(G,1), both(~ismember(both,ind_interf))), ~ismember(1:size(G,1), both(~ismember(both,ind_interf))));
                    pert_BoolData = BoolData(~ismember(1:size(G,1), both(~ismember(both,ind_interf))),:);

                    n=size(ind_interf,1); % perturbed nodes
                    k = dec2bin(0:(2.^n)-1)-'0'; % possible combinations

                    temp_res = zeros(size(k,1),4);
                    temp_name = cell(size(k,1),1);
		    temp_flipped = cell(size(k,1),1);

                    for l=1:size(k,1)
                        p_BoolData = pert_BoolData;
                        p_BoolData(p_BoolData(:,initial_state) ==2,initial_state) = k(l,:); % assign combinations of values
                        % Simulation
                        [Att_states, count] =boolSim_pert(pert_G,p_BoolData(:, initial_state),1000,ind_pert);
                        %Compare with unperturbed attractor
                        expected_Att_states = unperturbed_Att_states;
                        expected_Att_states(ind_pert) = ~expected_Att_states(ind_pert);
                        comparable_Att_states = Att_states(GRN_ind,:);
                        flipped_genes = sum(~(comparable_Att_states == expected_Att_states));
                        temp_res(l,:) = [sum(comparable_Att_states==expected_Att_states) / size(Att_states, 1)*100, ... 
                            sum(comparable_Att_states==expected_Att_states), count, flipped_genes]; 
			temp_flipped{l} = strjoin(GNames(find(comparable_Att_states ~= expected_Att_states)),';');

                        interf_name = strcat(GNames(ind_interf),'_',num2str(k(l,1)));
                        pert_name = strcat(GNames(ind_pert),'_',num2str(1*~unp_BoolData(ind_pert,initial_state)));
                        temp_name{l,:} = char(strjoin({interf_name{:},pert_name{:}},';'));
                    end

                end

                results = [results;temp_res];
                perturbations = [perturbations;cellstr(temp_name)];
		M_flipped = [M_flipped;cellstr(temp_flipped)];

            end
        end
    end       
end

%% Export results
fName = strcat(subfix,'-perturbation_4_interfaceNode.txt');  
fid = fopen(fName,'w');            
if fid ~= -1
    fprintf(fid,'%s\t%s\t%s\t%s\t%s\t%s\n','Perturbation', 'Match_fraction_initial', ...
        'Match_number_initial','Simulation_count','Flipped_genes','Names');
    for j=1:size(results,1)
        fprintf(fid,'%s\t',perturbations{j});       
        fprintf(fid,'%3.2f\t%d\t%d\t%d\t',results(j,1), results(j,2), results(j,3), results(j,4));      
	fprintf(fid,'%s\n',M_flipped{j});
    end
    % unperturbed
    fprintf(fid,'%s\t','Unperturbed');
    fprintf(fid,'%3.2f\t%d\t%d\t%d\t',m(1,1), m(1,2), m(1,3),m(1,4));     
    fprintf(fid,'%s\n',m_flipped);
    %fprintf(fid,'\n');
    fclose(fid);     
end

delete(MyPool);

